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//This is an implementation in C of a preferred embodiment 
//of the invention for the calculation of the Fourier transform 



#include <stdio.h> 
#include <math.h> 
#defme GEN 1000 
#define N 50 
#defineCORl 
#defmeEPS l.e-12 



//math.h is required only for the calculation of 0,+ 1 1 Fourier coeffs 
//the maximum allowed generation 
//output bandwidth 
//irrelevant parameter 

//tolerance, where contributions to Fourier coeffs of norm less than 
//EPS are regarded as negligible 



#define PI 3.141592653589793 



#defme NVAN 1 



//terminate the recursion when NVAN consecutive 
//generations have contributed negligibly 



#define MING 15 //calculate at least MING generations 



#define SMODE 1 



//SMODE=J includes the calculation of 0,+l,-l Fourier coeffs 
//SMODE=0 is faster and does not require math.h 



struct complex { //a complex number 

double x; 
double y; 

}; 

struct trip{ //convenient to have a real triple 

double alp; 
double bet; 
double gam; 

}; 

struct edge { //an arrow-structure 

int a,b,c,d; 
double e,alp,bet,gam; 

}; 

int count=0; // a running count of the total n umber of samples of input data 

struct complex fourier[2*N+l]; //the complex Fourier coeffs stored in the order 

//n=-N t ...,-l,0,l..-.,N, so fourier(n) is (N+n)th coeff 
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double abar,bbar,cbar; //coeffs in trig function abar*cos theta + bbar*sin theta + cbar 

//which agrees with input function for theta = -PI, -PI/2, 0 

void main(void){ 

void genertop(struct edge *,int,int,double) ; //recursive routine for top of circle 
void generbot(struct edge *,int,int,double); //recursive routine for bottom of circle 

void normalout(void); //output normalization routine 

void normalin(void); //input normalization routine 

double fbar(int,irtt); //fbar(p.q) returns f(theta) - (abar*cos theta + bbar*sin theta + cbar) 
//where theta is the argument of the complex number (p-iq)/(p+iq) 

struct edge doefl]; //initial edge for recursive calls 

int i; 

normalin(); //this calculates the normalization parameters abar, bbar, cbar 

doe->a=doe->d=l; //initialize the matrix of the doe to the identity 

doe->b=doe->c=0; 

doe->a!p=doe->bet=doe->gam=doe->e=0; //initialize wavelet and lagged trig coeffs 

genertop(doe,l,0,fbar(-l,!)); //begin recursion on top of the circle 

generbot(doe,l ,0,0.); //begin recursion on bottom of the circle 

normaloutQ; //this divides each output Fourier coeff by an overall factor PI 

//which was suppressed during the recursion and determines 
//negative from positive Fourier coeffs using reality of input 

printf("with EPS=%e, NVAN=%d, COR=%d, and MING=%d, the count is: %d\n", 
EPS,NVAN,COR,MING,count); 

for(i=0;i<=N;i++) 

printf("c%3d = %e+(i)%e\n",i,fourier[N+i].x,fourier[N+i].y); 



} 



void genertop(struct edge *p,int g.int envy .double fc){ //recursion for top of circle 

//input pointer to a current edge p of generation g, where there 
//have so far been envy consecutive generations of negligible 
//contributions, and fc is the current normalized input sample 
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edge twofer[2]; 



// the descendant edges of p 



struct complex temp; 



void genertop(struct edge * int,int,double); //recursive function 



struct complex tgener(struct edge *,struct edge * double); 



//tgener(p,twofer,f) stuffs twofer with the descendants of p 
//using the updated normalized input sample f returning the 
//next two normalized samples as its real and imaginary parts 



int charles(struct edge *); //charles(p) updates the Fourier coeffs with the 

//contribution from edge p returning 1 if they 
//are negligible, and returning 0 otherwise 

void tprune(struct edge *); //this implements the processing of a terminal edge 

if(envy>=NVAN&& g>=MING){ //if appropriate, then terminate cascade 
tprune(p); 
return; 
} 

temp=tgener(p,twofer,fc); //otherwise, generate descendants of p 

if(charles(twofer)*charles(twofer+l)) //if both descendants contribute negligibly, 

envy++; //then increment envy 



else 



envy=0; 



//otherwise reset envy to zero 



if(g>=GEN){ 



//problem: non-convergence after GEN generations 



printf{"fall through with gen=%d for a,b,c,d=%d,%d,%d > %d and e=%e\n", 
g,p->a,p->b,p->c,p->d,p->e); 

return; 



else{ 



//otherwise, continue recursion with descendants 



genertop(twofer,g+ 1 ,en vy.temp.x); 
genertop(twofer+ 1 ,g+ 1 ,en vy.temp.y); 
} 



void generbot{struct edge *p,int g.int envy .double fc){ 



//recursion for bottom of circle 
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struct edge twofer[2]; //this is entirely analogous to genertop 

struct complex temp; 

void generbot(struct edge *,int,int,doubIe); 

struct complex bgener(struct edge *,struct edge * double); 

int charles(struct edge *); 

void bprune(struct edge *); 

if(envy>=NV AN && g >= MING) { 
bprune(p); 
return; 
} 

temp=bgener(p,twofer,fc); 

if(charles(twofer)*charles(twofer+l)) 
envy++; 

else 

envy=0; 
if(g>=GEN){ 

printf("fall through with gen=%d for a,b,c,d=%d,%d,%d,%d and e=%e\n", 
g,p->a,p->b,p->c,p->d,p->e); 

return; 
} 

else{ 

generbot(twofer,g+l,envy,temp.x); 

generbot(twofer+l,g+l,envy,temp.y); 

} 



//charles(p) updates the Fourier coeffs using p 
//and returns a control parameter flag=l if 
//contribution is negligible and flag=0 otherwise 



it charles(struct edge * p){ 
int n,fiag; 

struct complex maJcecpx(int,int); //input p,q returns the complex number (p-iq)/(p+iq) 
struct complex mult(struct complex.struct complex); //complex multiplication 



int a,b,c,d; 
double e; 



struct complex pzl,pzr,pzf,pzb; 
struct complex pim,prm,pfm,pbm 



//local variables for matrix values 
//local variable for wavelet coeff 

//l,r,f,b labels left,right,front,back vertices 
//of Farey tesselation near the edge p, and 
//in a loop below on index n, pzx=(pxm) A n 
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//for x=l,r,f,b 



struct complex tmp,mul,cn; 



sixupc(int,int,int,int,double); 



//this updates the Fourier coeffs 0,+l,-l 
//provided SMODE=l 



a=p->a; 
b=p->b; 
c=p->c; 
d=p->d; 
e=p->e; 



//stuff local values 



pzr.x=b-d; 



//set-up for loop 



Fourier coeffs 



pzr.y=c-a; 

prm=makecpx(b-d,a-c) ; 
pzr=mult(pzr,pzr); 



pzl.x=b+d; 
pzl.y=-c-a; 

plm=makecpx(b+d,a+c); 
pzl=mult(pzl,pzl); 



pzf.x=d; 
pzf.y=-c; 

pfm=makecpx(d,c); 
pzf=mult(pzf,pzf); 



pzb.x=b; 
pzb.y=-a; 

pbm=makecpx(b,a); 
pzb=mult(pzb,pzb); 

flag=l; //initialize flag 

tmp.x=0; 

if(SMODE) 

flag=sixupc(a,b,c,d,e); //update 0,+ 1 ,-1 modes and re-set flag to 0 if any 
//of these contributions is non-negligible 

for(n=2;n<=N;n++){ //loop over positive Fourier coeffs n>l 

tmp.y=e/(n-n*n*n); 

pzr=mult(pzr,prm); //recursively calculate powers 

pzl=mult(pzl,pim); 
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pzf=mult(i.. .,pfm); 
pzb=mult(pzb,pbm); 



cn.x=-pzr.x+2*pzf.x+2*pzb.x-pzi.x; 
cn.y=-pzr.y+2*pzf.y+2*pzb.y-pzl.y; 

mul=mult(cn,tmp); 

fourier[n+N].x+=mul.x; //update nth Fourier coeff 
fourier[n+N].y+=mul.y; 

if(flag=l && mul.x*mul.x+mul.y*mul.y>=EPS) //if contribution to jith Fourier 
fiag=0; //coeff is non-negligible, 

//then set flag to zero 

} 

return(flag); 



int sixupc(int a M b.int c.int d.double e) { 



//this updates the 0.+1.-1 Fourier coeffs 
//using the input matrix entries a,b,c,d 
//and wavelet coeff e 



double tp.tm.tr.tl; //p.m.r.l labels the terminal,initial,right,!eft vertices of Farey 

//tesselation near edge, and tx denotes corresponding angle 

double dtemp; 



//tb=l on the top and tb=-l on the bottom of the circle 



//this is the local version of the flag in Charles and is the 
//value returned to charles by sixupc 



struct complex makecpx(int,int), temp; 
fl=l; //set the flag to 1 



tb=l.; 
if(b+c<0) 

tb=-l.; 

//calculate the various angles 

temp=makecpx(d,-c); 
tp=acos(tb*temp.x); 

temp=makecpx(b,-a); 



//default to the top of the circle 
//if on the bottom of the circle 
//then re-set tb=-l. 
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tm=acos(tb*temp.x.,, 

temp=makecpx(b+d,-(a+c)); 
tl=acos(tb*temp.x); 

temp=makecpx(b-d,c-a); 
tr=acos(tb*temp.x); 

dtemp=e*{ //contribution to 0th Fourier coeff 

2.*tp*(c*c+d*d)+2.*tm*(a*a+b*b) 
-tl*((a+c)*(a+c)+(b+d)*(b+d)) 
-tr*((a-c)*(a-c)+(b-d)*(b-d)) 

); 

if(dtemp*dtemp>=EPS) //if this contribution is non-negligible, 
fl=0; //then re-set the flag 

fourier[N].x+=dtemp; //up-date the 0th Fourier coeff 

temp.x=e*( //contribution to real part of 1st 

tp*(c*c-d*d)+tm*(a*a-b*b) //Fourier coeff 

+tr*((b-d)*(b-d)-(a-c)*(a-c))/2. 
+tl*((b+d)*(b+d)-(a+c)*(a+c))/2. 
); 

fourier[N+l].x+=temp.x; //update real part of 1st Fourier coeff 



temp.y=e*( //contribution to imag part of 1st 

2*(tp*c*d+tm*a*b) //Fourier coeff 

-tr*(b-d)*(a-c) 
-fl*(b+d)*(a+c) 
); 

fourier[N+l].y+=temp.y; 

if(temp.x*temp.x+temp.y*temp.y>=EPS) 
fl=0; 

retum(fl); 

} 



; complex mult(struct complex u,struct complex v) { //complex multiplication 

struct complex temp; 



//update imag part of 1st Fourier coeff 

//if this contribution is non-negligible, 
//then re-set the flag to zero 
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temp.x=u.x*v.x-u.) .y; 

temp.y=u.x*v.y+v.x*u.y; 

return(temp); 



struct complex makecpx(int p,int q){ //input p,q returns the complex number (p+iq)/(p+iq) 

struct complex temp; 
double den; 

den=p*p+q*q; 

temp.x=(p*p-q*q)/den; 
temp.y=-2*p*q/den; 

return(temp); 



void normalout(void) { //output normalization 

double fz=0.,fi=0.,fo=0.; 
double alp,bet,gam; 
int n; 

for(n=0;n<=N;n++) { 

fourier[N+n] .x=fourier[N+n].x/PI 
fourier[N+n].y=fourier[N+n].y/PI 
fourier[N-n].x=fourier[N+n].x; 
fourier[N-n].y=-fourier[N+n].y; 
) 



; //include factor PI suppressed before 

//calculate negative from positve 
//Fourier coeffs 



void normalin(void){ 

double fz=0.,fi=0.,fo=0.; 



//input normalization 

//values of input function for angles -PI, -PI/2, 0 



for(i=0;i<=2*N;i++) 

fourier[i].y=fourier[i].x=0; 



//f(p,q) returns the un-normalized input sample 
//at the angle corresponding to the point on the 
//circle given by the complex number (p-iq)/(p+iq) 
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fz=f(0,l); 
fi=f(l,0); 
fo=f(l,l); 

abar=(fi-fz)/2; //abar*cos theta + bbar*sin theta + cbar takes the same values 

cbar=(fi+fz)/2; //as f for the angles theta=-PI, -PI/2, 0 

bbar=cbar-fo; 



fourier[N].x=cbar*PI; 
fourier[N].y=0.; 



//stuff the 0th Fourier coeff scaled by PI to 
//recover it after re-scaling by I/PI in normalout 



fourier[N+ 1 ] .x=abar*PI/2; 
fourier[N+ 1 ] . y=-bbar*PI/2; 



//stuff the 1st Fourier coeff scaled by PI to 
//recover it after re-scaling by 1/PIin normalout 



double fbar(int p, int q){ 

struct complex temp; 
struct complex makecpx(int.int); 
double f (int,int); 
temp=makecpx(p,q); 

retum(f(p,q) - (abar*temp.x + bbar*temp.y+cbar)); 

} 



//input p,q to fbar produces the nornalized value 
//of the input data at the point (p-iq)/(p+iq) 



double f(int p,int q){ 

double fsinm(int,int,int); 
return(fsinm(p,q,6)); 

} 



//f is the input data, where f(p,q) is the 
//value taken at (p-iq)/(p+iq) 



double fsinm(int p.int q, int m){ //for example, the input function is 

1Rt n > //here taken to be sin(6* theta) 

struct complex zeta,meta; 

struct complex makecpx(int,int),mult(struct compiex,struct complex); 

zeta=makecpx(p,q); 

meta.y=0; 

for(n=l;n<=m;n++) 

meta=mult(zeta,meta); 



return(meta.y); 



struct complex tgener(struct edge *pc,struct edge *pn,double fc){ 

//tgener(p, twofer.fc) stuffs twofer with the descendants 
//of p for the top of the circle, where fc is the current 
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//normalized input sample, and tgener returns a complex nuiiioer. 
//the real and imag pans being the normalized input values 
//required for the least-squares regularization 



double fbar(int,int); 

int a,fa,c,d,bpd,apc; 

double e,alp,bet,gam; 

double axl,ayl,ax2,ay2; 
double bxl,byl,bx2,by2; 
double cxl,cyl,cx2,cy2; 

double sigl,taul,sig2,tau2,chi; 

double eu,et; 

double fn,fnp; 



//fbar(p,q) returns the normalized input values 

//local variables for matrix entries, bpd=b+d, apc=a+c 

//local variables for wavelet and lagged trig coeffs 

//index 1,2 refers to first.second (counter-clockwise) 
//descendants, prefix a.b.c refers to alp.bet.gam 
//update procedures below define indices x,y 

//quantities used in the least-squares regularization 

//wavelet coeffs of first and second descendant edges 

//normalized input samples 



double vlc,vlu,vlt,v2c,v2u,v2t; 

//coeffs in the inhomogeneous linear expression vic-fviu*eu+vit*et, 1=1,2, 
//of ongoing approximation to normalized input at the respective sample points 
//((2b+d)+i(2a-fc))/((2b+d)-i(2a+c)),«b+2d)+i{a+2c))/((b+2d)-i(a+2c)) 



struct complex temp; 
struct edge *pnp; 



pnp=pn; 
pnp++; 



//initialize pointers 



a=pc->a; 
b=pc->b; 

d=pc->d; 
apc=a+c; 
bpd=b+d; 



//initialize matrix entries 



pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 



//stuff matrix of first descendant 



pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 



//stuff matrix of second descendant 
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alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



//initialize lagged trig coeffs 



e=pc->e; 



//initialize wavelet coeff 



temp.x=fn=fbar(-(b+bpd),a+apc); //store new normalized input samples as the 
temp.y=fnp=fbar(-(d+bpd),c+apc); //real and imag parts of temp to return 



axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 
ayl=2*(b*b-a*a); 

bx l=bet+(c*d-a*d-b*c-a*b)*4*e; 
byl=4*a*b; 

cx 1 =gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; 
cyl=-2*(a*a+b*b); 

ax2=aip+4*(d*d-c*c)*e; 
ay2=2*(c*c-d*d); 

bx2=bet+8*e*c*d; 
by2=-4*c*d; 

cx2=gam-4*e*(c*c+d*d); 
cy2=2*(c*c+d*d); 

chi=2*e+( //chi=eu-et 

((a+c)*(a+c)+(b+d)*{b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

} 

else i //prepare for stuffing alp,bet,gam in case c>=a 



axl=alp+4*e*(a*a-b*b); 
ayl=2*(b*b-a*a); 

bxl=bet-8*e*a*b; 
byi=4*a*b; 

cxl=gam+4*e*(a*a+b*b); 



count+=2; 



//update count of total sample points 



if(a>c){ 



//prepare for stuffing alp.bet.gam in case a>c 
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cyi=-2*(a-a+b*b); 

ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c+c*d)*4*e; 
by2=-4*c*d; 

cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; 
cy2=2*(c*c+d*d); 



chi=-2*e+( //chi=eu-et 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

} 

//in any case, calculate the coeffs vic.viu.vit, for i=l,2 
vlc=-fn+cxl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl 
+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+-bpd)+(a+apc)*(a+apc)); 

vlt=cyl+ 

(<(b+bpd)*<b+bpd)-(a+apc)*(a+apc))*ayl 
+2*byl*(b+bpd)*(a+apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

v 1 u=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

v2c=-fnp+cx2+ 

<((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2 
+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2 
+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

sigl=vlc+chi*vlu; 
taul=vlt+vlu; 

sig2=v2c+chi*v2u; 
tau2=v2t+v2u; 

et=-(sigl*taul+sig2*tau2)/(taul*taul+tau2*tau2); //calculate et and eu 
eu=chi+et; 
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pn->e=eu; 
pnp->e=et; 



//stuff: 



wavelet coeffs 



pn->alp=ax 1+ayl *et; 
pn->bet=bxl+byl*et; 
pn->gam=cxl+cy 1 *et; 



//stuff first lagged trig coeffs 



pnp->a!p=ax2+ay2*eu; 
pnp->bet=bx2+by2*eu; 
pnp->gam=cx2+cy2*eu; 



//stuff second lagged trig coeffs 



return (temp); 



struct complex bgener(struct edge *pc,struct edge *pn,double fc){ 



//bgener(p,twofer,fc) stuffs twofer with the descendants 
//of p for the bottom of the circle, where fc is the current 
//norma!i2ed input sample and returns a complex number, 
//the real and imag parts being the normalized input values 
//required for the least-squares regularization 

//this is entirely analogous to tgener, and only the differences 
//will be noted here 



double fbar(int,int); 

int a,b,c,d,bpd,apc; 
double e,alp,bet,gam; 

double axi,ayl,ax2,ay2; 
double bxl,byl,bx2,by2; 
double cxl,cyl,cx2,cy2; 

double sig 1 ,tau 1 ,sig2,tau2,chi; 
double vlc,vlu,vlt,v2c,v2u,v2t; 

double eu,et,fn,fnp; 

struct complex temp; 

struct edge *pnp; 

pnp=pn; 
pnp++; 
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a=pc->d; //differs from tgener in that d,-c,-b,a replaces a.b.c.u 

b=-(pc->c); 

c=-(pc->b); 

d=pc->a; 

apc=a+c; 
bpd=b+d; 

pn->a=bpd; //differs from tgener in that d,-c,-b,a replaces a,b,c,d 

pn->b=-apc; 

pn->c=-b; 

pn->d=a; 

pnp->a=d; //differs from tgener in that d,-c,-b,a replaces a,b,c,d 

pnp->b=-c; 

pnp->c=-bpd; 

pnp->d=apc; 

aip=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



temp.x=fn=fbar(a+apc,b+bpd); //differs from tgener in that 
temp.y=fnp=fbar(c+apc,d+bpd); //(-q,p) replaces (p,q) 



axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 
ayl=2*(b*b-a*a); 

bxl=bet+(c*d-a*d-b*c-a*b)*4*e; 
byl=4*a*b; 

cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; 
cyl=-2*(a*a+b*b); 

ax2=alp+4*(d*d-c*c)*e; 
ay2=2*(c*c-d*d); 

bx2=bet+8*e*c*d; 
by2=-4*c*d; 

cx2=gam-4*e*(c*c+d*d); 
cy2=2*(c*c+d*d); 
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chi=2*e+( 

«a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

J 



else{ 

axl=alp+4*e*(a*a-b*b); 
ayl=2*(b*b-a*a); 

bxl=bet-8*e*a*b; 
byl=4*a*b; 

cxl=gam+4*e*(a*a+b*b); 
' cyl=-2*(a*a+b*b); 

ax2=alp+{a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c+c*d)*4*e; 
by2=-4*c*d; 

cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; 
cy2=2*(c*c+d*d); 



chi=-2*e+( 

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam) 

-alp*((b+d)*(b+d)-(a+c)*(a+c)) 

-bet*2*(a+c)*(b+d) 

)/4.; 

} 

vlc=-fn+cxl+ 

(((b+bpd)*(b+bpd)-(a+apc)*{a+apc))*axl 
+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

vlt=cyl+ 

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl 
+2*byl*(b+bpd)*(a+apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 
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v2c=-fnp+cx2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c-t-apc))*ax2 
+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+ 

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2 
+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc»; 

v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

sigl=vlc+chi*vlu; 
taul=vlt+vlu; 

sig2=v2c+chi*v2u; 
tau2=v2t+v2u; 

et=-(sigl*taul+sig2*tau2)/(taul*taul+tau2*tau2); 
eu=chi+et; 

pn->e=eu; 
pnp->e=et; 

pn->alp=axl+ayl*et;; 
pn->bet=bx 1+by 1 *et; 
pn->gam=cxl+cyl*et; 

pnp->alp=ax2+ay2*eu; 

pnp->bet=bx2+by2*eu; 
pnp->gam=cx2+cy2*eu; 

return (temp); 



void tprune(struct edge *p) { //this routine implements the processing of terminal 

//edge p in the cascade for the top of the circle 

double alp, bet, gam.e; //local copies of lagged trig coeffs and wavelet coeff 

double alpl,alp2,betl,bet2,gaml,gam2; //trig coeffs of the trigonometric function sigma at 
//the first and second (counter-clockwise) samples 

struct trip mob(struct complex,struct complex, struct complex, 
double,double,doubie); 
//mob(zl,z2,z3,fl,f2,f3) returns the triple with entries alp.bet.gam 
//where alp*cos theta + bet*sin theta + gam takes the respective values 
// fl,f2,f3 at the points zl,z2,z3 on the circle given as complex numbers 
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double fbar(int,int); //returns normalized input values 

void fixup(struct trip,int,int,int.int,int); //adjusts the output Fourier coeffs using sigma-tau 

struct complex makecpx(int.int); //input p,q returns (p-iq)/(p+iq) 

void sixupp(struct trip,int,int,int,int,int); //adjusts the 0.+1.-1 Fourier coeffs provided 
//SMODE=l 

int a,b,c,d; //local variables for the matrix entries 

struct trip tripl,trip2; //first and second trig triples 

a=p->a; //stuff local variables 

b=p->b; 

c=p->c; 

d=p->d; 

e=p->e; 

alp=p-:>alp; 

bet=p->bet; 

gam=p->gam; 

if (a>c){ //final update of trig fields in case a>c 

alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*<b+d)-c*(c-a)); 

beti=bet+e*4*(c*(d-b)-a*(d+b)); 

gaml=gam+e*2*(a*(a-K:)+c*(a-c)+b*(b+d)-d*(d-b)); 

aip2=alp+e*4*(d-c)*(d+c); 
bet2=bet+e*8*c*d; 
gam2=gam-e*4*(c*c+d*d); 
) 

elsef //final update of trig fields in case o=a 

alpl=alp+e*4*(a-b)*(a+b); 
betl=bet-e*8*a*b; 
gam 1 =gam+e*4*(a*a+b*b); 



alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); 
bet2=bet+e*4*(a*(d-b)+c*(b+d)); 
gam2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d)); 
) 

tripl=mob( //calculate sigma for first descendant 

makecpx(-(3 *b+d),3 *a+c), 
makecpx(-(2*b+d),2*a+c), 
makecpx(-(3*b+2 *d),3*a+2*c), 
fbar(-(3*b+d),3*a+c), 
fbar(-(2*b+d),2*a+c), 
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fbar(-(3*b-^*d),3*a+2*c) 

); 

trip2=mob( //calculate sigma for second descendant 

makecpx(-(2*b+3*d),2*a+3*c), 
makecpx(-{b+2*d),a+2*c), 
makecpx(-(b+3*d),a+3*c), 
fbar(-(2*b+3 *d),2 *a+3 *c), 
fbar(-(b+2*d),a+2*c), 
fbar(-(b+3*d),a+3*c) 
); 

count+=6; //update count of total sample points 

tripl.alp=(tripl.alp-alpl)/COR; //difference of first trig functions 

tripl.bet=(tripl.bet-betl)/COR; 

trip 1 .gam=(trip 1 .gam-g am 1 )/COR ; 

trip2.alp=(trip2.alp-alp2)/COR; //difference of second trig functions 

trip2.bet=(trip2.bet-bet2)/COR; 

trip2.gara=(trip2.gam-gam2)/COR; 

fixup(tripl,b+d,a+c,b,a,l); //adjust Fourier coeffs with first trig function difference 

fixup(trip2,d,c,b+d,a+c,l); //adjust Fourier coeffs with second trig function difference 

if(SMODE){ //if SMODE=l, then adjust O.H-1,-1 Fourier coeffs as well 

sixupp(trip 1 ,b+d,a+c,b,a, 1 ); 
sixupp(trip2,d,c,b+d,a+c, 1); 
} 



void bprune(struct edge *p){ //this routine implements the processing of terminal 

//edge p in the cascade for the bottom of the circle 

//this is entirely analogous to tgener, and only the 
//differences will be noted here 

double alp, bet, gam.e; 
double alp 1 ,alp2 ,be 1 1 ,bet2,gam 1 ,gam2; 
struct trip mob(struct complex.struct complex, struct complex, 
double,double,double); 

double fbar(int,int); 

void fixup(struct trip, int,int,int,int,int); 
struct complex makecpx(int,int); 
void sixupp(struct trip,int,int,int,int,int); 
int a,b,c,d; 

struct trip trip l,trip2; 
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a-p->d; //differs from tprune in that d,-c,-b,a replaces a,b,c,d 

b=-(p->c); 

c=-(p->b); 

d=p->a; 

e=p->e; 

alp=p->aip; 

bet=p->bet; 

gam=p->gam; 

if(a>c){ 

alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); 
bet 1 =bet+e*4*(c*(d-b)-a*(d+b)); 
gami=gam+e*2*(a*(a+c)+c*(a-c)+b*(b+d)-d*(d-b)); 

alp2=alp+e*4*(d-c)*(d+c); 
bet2=bet-i-e*8*c*d; 
gam2=gam-e*4*(c*c+d*d); 
} 

else{ 

alpl=aJp+e*4*(a-b)*(a+b); 

betl=bet-e*8*a*b; 

gaml=gam+e*4*(a*a+b*b); 

alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); 
bet2=bet+e*4*(a*(d-b)+c*(b+d)); 
gam2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d)); 
} 



tripl-mob( //arguments of mob differ from those in tprune 

makecpx(3*a+c,3*b+d), 
makecpx(2*a+c,2*b+d), 
makecpx(3 *a+2*c,3 *b+2*d), 
fbar(3*a+c,3*b+d), 
fbar<2*a+c,2*b+d), 
fbar(3 *a+2*c,3*b+2*d) 

); 

trip2=mob( //arguments of mob differ from those in tprune 

makecpx(2*a+3*c,2*b+3*d), 
makecpx(a+2*c,b+2*d), 
makecpx(a+3*c,b+3*d), 
fbar(2*a+3*c,2*b+3*d), 
fbar(a+2*c,b+2*d), 
fbar(a+3*c,b+3*d) 
); 

count+=6; 
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tripl.alp=(tripl.alp+alpl)/COR; //differs from tprune in that alp.bet 
tripl.bet=(tripl.bet+betl)/COR; //is replaced by -alp,-bet 
trip 1 .gam=(tripl .gam-gam 1 )/COR; 

trip2.alp=(trip2.alp+aip2)/COR; //differs from tprune in that alp,bet 
trip2.bet=(trip2.bet+bet2)/COR; //is replaced by -alp,-bet 
trip2.gam=(trip2.gam-gam2)/COR; 

fixup(tripl,b+d,a+c,b,a,-l); //last argument of fixup differs from tprune 

fixup(trip2,d,c,b+d,a+c,-l);Q //last argument of fixup differs from tprune 

if(SMODE){ //if SMODE=l, adjust the 0.+1.-1 Fourier coeffs 

sixupp(tripl ,b+d,a+c,b,a,-l); //last argument of sixupp differs from tprune 
sixupp(trip2,d,c,b+d,a+c,-l); //last argument of sixupp differs from tprune 



struct trip mob(struct complex zl,struct complex z2,struct complex z3, 
double fl,double O.double B){ 
//mob{zl,z2,z3,fl,f2,f3) returns the triple with entries alp,bet,gam 
//where afp*cos theta + bet*sin theta + gam takes the respective values 
// fl,f2,f3 at the points zl,z2,z3 on the circle given as complex numbers 

double det,aa,bb,cc,dd; 
struct trip temp; 

aa={zl.x-z3.x); 
bb=(zl.y-z3.y); 
cc=(z2.x-z3.x); 
dd=(z2.y-z3.y); 

det=aa*dd-bb*cc;//general principles guarantee that det is non-zero 

aa=aa/det; 
bb=bb/det; 
cc=cc/det; 
dd=dd/det; 

temp.aIp=dd*(fl-f3)-bb*{f2-0); 
temp.bet=aa*(f2-f3)-cc*(fl-f3); 
temp.gam=0-z3.x*temp.alp-z3.y*temp.bet; 

return(temp); 
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void fixup(struct trip delt.int dm.int cm,int dp.tnt cp.int tb){ 

//fixup delt,dm,cm,dp,cp,tb) adjusts the Fourier coeffs by adding the Fourier 
//coeffs of delt=T A C-T with support truncated to be the interval with endpoints 
//(dm+i *cm)/(dm-i *cm), (dp+i *dp)/(dp-i *cp) 

/Aet tm, tp be the corresponding angles 

// tb=+l on the top, and tb=-l on the bottom of the circle 

int n; 

double cpp,cmp,cpm,cmm; 
//in loop below on n, cpx=[cos(l+n)tx]/{n+l), cmx=[cos(l-n)tx]/(l-n), for x=p,m 

double spm,smm, spp,smp; 
//in loop below .on n, spx=[sin(l+n)tx]/(l+n). smx=[sin(l-n)tx]/(l-n), for x=p,m 

double crm,crp,srm,srp; 
//in loop below on n, crm=[cos n*tx]/n, srm= [sin n*tx]/n, for x=p,m 

double alp,bet,gam; //local variables for entries of delt 

double ctp, ctm, cnp.cnm; //in loop below on n, ctx=cos tx, cnx=cos n*tx, for x=p,m 
double stm,snm,stp,snp; //in loop below on n, stx=sin tx, snx=sin n*tx, for x=p,m 

double cbufp,sbufp,cbufm,sbufm; 
struct complex makecpx(int,inc),temp; 



alp=delt.alp; //stuff local variables 

bet=delt.bet; 

gam=delt.gam; 

temp=makecpx(dm,-cm); 

ctm=tb*temp.x; 

stm=tb*temp.y; 

temp=makecpx(dp,-cp); 

ctp=tb*temp.x; 

stp=tb*temp.y; 

cnp=ctp*ctp-stp*stp; 
snp=2*stp*ctp; 

cnm=ctm*ctm-stm*stm; 
snm=2*stm*ctm; 
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for(n=2;n<=N;n++){ //loop over Fourier coeffs 

cbufp=ctp*cnp-stp*snp; //update using formulas for cos 
sbufp=stp*cnp+ctp*snp; //or sin of a sum of angles 
cbufm=ctm*cnm-stm*snm; 
sbufm=stm*cnm+ctm*snm; 

cpp=cbufp/(l+n); //update for p 

cmp=(cbufp+2*stp*snp)/( 1 -n); 
spp=sbufp/(l+n); 
smp=(sbufp-2*ctp*snp)/(I -n); 

cpm=cbufm/(l+rt); //update for m 

cmm=(cbufm+2*stm*snm)/(l-n); 
spm=sbufm/(l+n); 
smm=(sbufm-2*ctm*snm)/{l-n); 

crp=2.*cnp/n; 
crm=2.*cnm/n; 
. srp=2.*snp/n; 
srm=2.*snm/n; 

cnp=cbufp; 
snp=sbufp; 
cnm=cbufm; 
snm=sbufm; 

fourier[N+n].x+=( 
+alp*(smp+spp-smm-spm) 
+bet*(-cpp-cmp+cpm+cmm) 
+gam*(srp-srm) 
)/4.; 

fourier[N+n].y+=( //adjust imag part of nth Fourier coeff 

+alp*(cpp-cmp-cpm+cmm) 

+bet*(spp-smp-spm+smm) 

+gam*(crp-crm) 

)/4.; 



//update for next iteration 



//adjust real part of nth Fourier coeff 



void sixupp(struct trip triple.int dm.int cm.int dp.int cp.int tb)( 

//sixup(delt,dm,cm,dp,cp,tb) adjusts the 0,+l,-l Fourier coeffs by adding the 0.+1.-1 Fourier 
//coeffs of de!t=T A C-T with support truncated to be the interval with endpoints 
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//(dm+i *cm)/Cdm-i *cm), (dp+i -*dp)/(dp-i *cp) 

//let tm, tp be the corresponding angles 

// tb=+l on the top, and tb=-l on the bottom of the circle 



double aip,bet,gam; //local variables for entries of delt 

double com,sim,cop,sip,tm,tp; //cox=cos tx, six=sin tx, for x=p,m 

struct complex temp; 

struct complex makecpx(int,int); 



temp=makecpx(dm,-cm); //calculate com.sim.tm 

com=tb*temp.x; 

sim=tb*temp.y; 

tm=acos(com); 

temp=makecpx(dp,-cp); //calculate cop,sip,tp 

cop=tb*temp.x; 

sip=tb*temp.y; 

tp=acos(cop); 



alp=tb*triple.alp; //adjust alp.bet for the bottom of the circle 

bet=tb*triple.bet; 

gam=triple.gam; 

fourierfN] .x+=( //update 0th Fourier coeff 

gam*(tp-tm)+alp*(sip-sim)-bet*(cop-com) 
)/2.; 

fourier[N+l].x+=( //update real part of 1st Fourier coeff 

aIp*(tp-tm)/2. 
+gam*(sip-sim) 

-bet*(cop*cop-sip*sip-com*com+sim*sim)/4. 

+alp*(cop*sip-com*sim)/2. 

)/2.; 

fourier[N+l].y+=( //update imag part of 1st Fourier coeff 

-bet*(tp-tm)/2. 
+gam*(cop-com) 

+alp*(cop*cop-stp*sip-com*com+-sim*sirn)/4. 

+bet*(cop*sip-com*sim)/2. 

)/2.; 
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//Here are the subroutines replacing tgener and bgener above in an 
//implementation which takes advantage of renormalization 

struct complex tgener(struct edge *pc,struct edge *pn,doubIe fc) 
{ 



double fbar(int,int); 
int a,b,c,d,bpd,apc; 
double e,alp,bet,gam; 
double axi,ayl,ax2,ay2; 
double bxl,byl,bx2,by2; 
double cxl,cyl,cx2,cy2; 
double sig 1 ,tau 1 ,sig2,tau2,chi; 
double vlc,vlu,vlt,v2c,v2u,v2t; 
double deriv(int,int,int,int,int,int); 
double dc,dn,dnp; 
double eu,et,fn,fnp; 
struct complex temp; 
struct edge *pnp; 
extern int count; 

pnp=pn; 
pnp++; 



a=pc->a; 
b=pc->b; 
c=pc->c; 
d=pc->d; 

apc=a+c; 
bpd=b+d; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 



49 



WO 00/44104 



PCT/US99/30584 



a!p=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 



temp.x=fn=fbar(-(b+bpd),a+apc); 
temp.y=fnp=fbar(-(d+bpd),c+apc); 
temp.n=0; 

count+=2; 

dc=deri v(a,b,c,d,- 1,1); 
dn=deri v(a,b,c,d,- 1 ,2); 
dnp=deriv(a,b,c,d,-2, 1); 

fc=fc/dc; 
fn=fn/dn; 
fnp=fnp/dnp; 



if(a>c){ 

axl=alp+4.*e; 

bxi=bet-4.*e; 

cxl=gam; 

ax2=alp+4 *e; 

bx2=bet; 

cx2=gam-4.*e; 

chi=2.*((fc-gam-bet)/4.+e); 

} 



else{ 



ax l=alp+4.*e; 

bxl=bet; 

cxl=gam+4.*e; 

ax2=alp+4.*e; 

bx2=bet+4.*e; 

cx2=gam; 

chi=2.*((fc-gam-bet)/4.-e); 
} 

v 1 c=-f n+cx 1 +(4.*bx 1 -3.*ax 1 )/5. ; 

vlt=-4./5.; 

viu=8./5.; 
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v2c=-fnp+cx2+(4. * ^2+3 . *ax2)/5 . ; 

v2u=4./5.; 

v2t=-8./5.; 

sigl=vlc+chi*vlu; 
' taul=vlt+vlu; 

sig2=v2c+chi*v2u; 
tau2=v2t+v2u; 

et=-(sigl*taul*(l)+sig2*tau2*(l))/(taul*taul*(l) 
+tau2*tau2*(l)); 

eu=chi+et; 

pn->e=eu; 
pnp->e=et; 

alp=axl-2.*et; 

bet=bxl; 

gam=cxl-2.*et; 

pn->alp=(2.*bet+alp+gam)/2.; 

pn->bet=(bet-alp+gam); 

pn->gam=(2.*bet-aIp+3.*gam)/2.; 

alp=ax2-2.*eu; 

bet=bx2; 

gam=cx2+2.*eu; 

pnp->alp=(-2.*bet+alp-gam)/2.; 

pnp->bet=alp+bet+gam; 

pnp->gam==(2.*bet+a!p+3.*garn)/2.; 



return (temp); 



double deriv(int a.int b.int c.int d,int p.int q){ 

double x.y.nep; 

x=p*p-q*q; 
y=-2*p*q; 
X=x/(p*p+q*q); 

y=y/(p*P+q*q); 

nep=-(a*a+b*b+c*c+d*d)+x*(a*a+b*b-c*c-d*d)-y*(2*a*c+2*b*d); 
return(-2./nep); 
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struct complex bgener(struct edge *pc,struct edge *pn,doub)e fc) 
f 

double fbar(int,int); 
int a,b,c,d,bpd,apc; 
double e,alp,bet,gam; 
double axl,ayl,ax2,ay2; 
double bxl,byl,bx2,by2; 
double cxl,cyl,cx2,cy2; 
double sigl,taul,sig2,tau2,chi; 
double vlc,vlu,vlt,v2c,v2u,v2t; 
double eu,et,fn,fnp; 
extern int count; 

double deriv(int,int,int,int,int,int); 
double dc.dn.dnp; 
struct complex temp; 
struct edge *pnp; 

pnp=pn; 
pnp++; 

a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 

apc=a+c; 
bpd=b+d; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 

pnp->a=d; 
pnp->b=-c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 

temp.x=fn=fbar(a+apc,b+bpd) ; 
temp.y=frtp=fbar(c+apc,d+bpd); 
temp.n=0; 
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count+=2; 

dc=deriv(d,-c,-b,a,l,l); 
dn=deriv(d,-c,-b,a,2,I); 
dnp=deriv(d,-c,-b,a,l,2); 

fc=fc/dc; 
fn=fn/dn; 
fnp=fnp/dnp; 

if(a>c){ 

axl=a!p+4.*e; 
bxl=bet-4.*e; 
cxl=gam; 

ax2=alp+4.*e; 

bx2=bet; 

cx2=gam-4.*e; 

chi=2.*((fc-gam-bet)/4.+e); 
} 



else{ 

axl=alp+4*e; 

bxl=bet; 

cxl=gam+4*e; 

ax2=alp+4*e; 
bx2=bet44*e; 
cx2=gam; 

chi=2*((fc-gam-bet)/4-e); 
} 

vlc=-fn+cxl+(4.*bxl-3.*axl)/5.; 

vlt=-4./5.; 

vlu=8./5.; 

v2c=-fnp+cx2+(4.*bx2+3.*ax2)/5.; 

v2u=4./5.; 

v2t=-8./5.; 

sigl=vlc+chi*vlu; 
taul=vlt+vlu; 

sig2=v2c+chi*v2u; 
tau2=v2t+v2u; 
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et=-(sig 1 *tau 1 *( 1 ) . _. g 2*tau2*( 1 ))/(tau 1 *tau 1 *(I ) 
+tau2*tau2*(l)); 

eu=chi+et; 
pn->e=eu; 
pnp->e=et; 

alp=axl-2.*et; 

bet=bxl; 

gam=cxl-2.*et; 

pn->alp=(2.*bet4-alp+gam)/2.; 

pn->bet=(bet-alp+gam); 

pn->gam=(2.*bet-alp+3.*gam)/2.; 

alp=ax2-2.*eu; 

bet=bx2; 

gam=cx2+2.*eu; 

pnp->alp=(-2. *bet+alp-gam)/2. ; 

pnp->bet=alp+bet+gam; 

pnp->gam=(2.*bet+alp+3.*gamy2.; 

return (temp); 
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//This is an implen .tation in C of a preferred embodiment 

//of the invention for the calculation of the inverse Fourier transform 



#include <stdio.h> 
#deftneSCAT50 
#defineN50 
#defineMINMODE0 
#define PI 3.141592653589793 
#define CHUNK 5000 



//SCAT=[exp sinh A (-l)(l/SCALE)J A 2 
//input bandwidth 
//minimum frequency of input 

//chunk of memory allocated as required 



struct complex { 



double x; 
double y; 



//a complex number 



//a complex arrow-structure 



struct edge{ 

int a,b,c,d; 
struct complex e,alp,bet,gam; 

}; 

struct pqxy{ //the basic data type of the output, where 

int p,q; //p,q corresponds to the complex number (p-iq)/(p+iq) 

double x,y; //at which the output function takes the complex value x+iy 
}; 



struct complex pfourier[2*N+l]; 
struct complex *fourier=pfourier+N; 

struct complex czer.cone.cmon; 
struct pqxy graf[CHUNK]; 

int outcount=0; 



//the input Fourier coeffs 
//indexed -N -1,0,1,...,N 

//the modified 0.+1.-1 Fourier coeffs 

//output spatial values in 

//clockwise-order around the circle 
//starting from the complex number -1 

//running tally of the total 
//number of output values 



voidmain(void){ 



void genertop(struct edge *,int); //recursive routine for top of circle 
void generbot(struct edge *,tnt); //recursive routine for bottom of circle 



void normalout(void); 
void normalin(void); 



//output normalization routine 
//input normalization routine 



struct complex cogi(struct edge *); //returns complex wavelet coeff of argument 
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struct complex ei,emu,emt,eu,et: //wavelet coeffs on I,U A {-1 },T A {-1 },U,T 

struct complex emtt,emtu,emut,emuu; //wavelet coeffs on T A {-2}, 

//T A {-1 }U A {-1 },U A {-1}T A {-1 },U A {-2} 

struct edge startfl]; //initial edge for each of 

//the various recursive calls 

int i; 

//begin input data 

for(i=0;i<=2*N;i++){ 

pfourier[i].x=0.; 
pfourier}i].y=0.; 
} 

fourier[5].x=.5; 
fourier[-5].x=.5; 

//end input data 

normaIin(); //this stuffs czer,cone,cmon with the 

//modified O.+l.-l Fourier coefficients 

graf[outcount].p=0; //initial normalized output value 

graf[outcount].q=I; 

graf[outcount].x=0.; 

graf[outcount].y=0.; 

outcount++; 



//calculate the wavelet coeffs for the various required edges of small generation 



start->a=start->d=l; 
start->c=0; 
start->b=-l; 
emt=cogi(start); 

start->a=start->d=l; 
start->c=0; 
start->b=l; 
et=cogi(start); 

start->a=start->d=l; 
start->c=-l; 
start->b=0; 
emu=cogi(start); 

start->a=start->d=l ; 
start->c=l; 
start->b=0; 
eu=cogi(start); 

start->a=start->d= 1 ; 
start->b=start->c=0; 
ei=cogi(start); 
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stan->a=start-xi= 1 ; 
start->b=0; 
start->c=-2; 
emuu=cogi{start); 

start->a=start->d=l; 
start->b=-2; 
start->c=0; 
emtt=cogi(start); 

start->a=I; 
start->b=start->c=-l ; 
start-x3=2; 
emut=cogi(start); 

start->a=2; 
start->d=l; 
start->b=start->c=-l ; 
emtu=cogi(start); 



//serially perform the recursions for the initial edges in counter-clockwise order 



start->a=start->d=l; //initializ 

start->c=l; 

start->b=0; 

start->e=eu; 

start->alp.x=2.*ei.x-2.*et.x; 

start->alp.y=2.*ei.y-2.*et.y; 

start->bet.x==2.*(emt.x-emu.x)-2.*ei.x; 

start->bet.y=2.*(emt.y-emu.y)-2.*ei.y; 

start->gam.x=2.*ei.x-2.*et.x; 

start->gam.y=2.*ei.y-2.*et.y; 

genertop(start, 1); //recursion for U 

start->a=start->d=l; //initializ. 

start->c=0; 

start->b=l; 

start->e=et; 

start->alp.x=2.*ei.x-2.*eu.x; 

start->alp.y=2.*ei.y-2.*eu.y; 

start->bet.x=2.*(emt.x-emu.x)+2.*ei.x; 

start->bet.y=2.*(emt.y-emu.y)+2.*ei.y; 

start->gam.x=-2.*ei.x+2.*eu,x; 

start->gam.y=-2.*ei.y+2.*eu.y; 

genertop(start,I); //recursion for T 

start->a=start->d=l; //initialize 

start->c=0; 

start->b=-2; 

start->e=emtt; 

start->alp.x=-(-2.*ei.x+2.*emu.x-4.*emt.x+2.*emut.x); 
start->alp.y=-(-2.*ei.y+2*emu.y-4.*emt.y+2.*emut.y); 
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start->bet.x=-(2. *e, .. . -2*erau.x+2.*emt.x) ; 
stan->bet.y=-(2.*ei.y-2.*emu.y+2.*emt.y); 
start->gam.x=2.*ei.x-2 *emu.x+4.*emt.x-2 *emut.x; 
start->gam.y=2.*ei.y-2.*emu.y+4.*emt.y-2.*emut.y; 



generbot(start,l); 



//recursion for T A { -2} 



start->a=l; 



//initialize 



start->d=2; 
start->b=start->c=- 1 ; 
start->e=emut; 

start->alp.x=-(-2.*eijc+2.*ernu.x+2.*emt.x); 

start->alp.y=-(-2.*ei.y+2.*emu.y+2.*emt.y); 

start->bet.x=-(2.*ei.x-2.*emu.x-6.*emt.x+4.*emtt.x); 

start->faet.y=-(2.*ei.y-2.*emu.y-6.*emt.y44.*emtt.y); 

start->gam.x=2.*ei.x-2.*erau.x-6.*emt.x+4.*emtt.x; 

start->gam.y=2.*ei.y-2.*emu.y-6.*emt.y+4.*emtt.y; 

generbot(start,l); //recursion for U A {-1 }T A {-1 } 



start->a=2; //initialize 

start->d=l; 

start->b=start->c=-l; 

start->e=emtu; 

start->alp.x=-(-2.*ei.x+2.*emu.x+2.*emt.x); 

start->alp.y=-(-2.*ei.y+2.*emu.y+2.*emt.y); 

start->bet.x=-(-2.*ei.x+6.*emu.x+2.*emt.x-4.*emuu.x); 

start->bet.y=-(-2.*ei.y+6.*emu.y+2.*emt.y-4.*emuu.y); 

start->gam.x=-2.*ei.x+6.*emu.x+2.*emt.x-4.*emuu.x; 

start->gam.y=-2.*ei.y+6.*emu.y+2.*emt.y-4.*emuu.y; 

generbot{start,l); //recursion forT A {-l }U A {-1 } 

start->a=start->d=l; //initialize 

start->b=0; 

start->c=-2; 

start->e=emuu; 

start->alp.x=-(-2.*ei.x-4*emu.x+2.*emt.x+2.*emtu.x); 

start->alp.y=-(-2.*ei.y-4.*emu.y+2.*emt.y+2.*emtu.y); 

start->bet.x=-(-2.*ei.x-2.*emu.x+2 *emt.x); 

start->bet.y=-(-2.*ei.y-2.*emu.y+2.*emt.y); 

start->gam.x=-2.*ei.x-4.*emu.x+2.*emt.x+2.*emtu.x; 

start->gam.y=-2.*ei.y-4.*emu.y+2.*emt.y+2.*emtu.y; 

generbot(starU); //recursion for U A {-2} 

//recursions complete, and graf is stuffed with normalized output data in counter-clockwise order 

normaloutO; //un-normalize by altering 



//each output value using the 
//modified O.+l.-l Fourier coeffs 



printf("cmon=%e +(i) %e\n czer=%e +(i) %e\n cone=%e +(i) %e\n\n", 
cmon.x,cmon.y,czer.x,czer.y,cone.x,cone.y); 



for (t=0;i<outcount;i++) 
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printf("(%^. jd) %e+(i)%e\n", 

graf[i].p,grafli].q,graf[i].x,grafli].y); 

} 



void genertopCstruct edge *p,int g){ //recursion for top of the circle 

struct edge twofer[2],*ped; //recursion generates twofer[2] fromp, where 



struct complex temp; 

struct complex makecpx(int.int); 

double ct,st; 

void genertop(struct edge *,int); 

void tgener(struct edge *,struct edge *); 



int subtend(struct edge *),pp,qq; 



tgener(p,twofer); 



if(subtend(twofer)){ 



:s of twofer are the descendants of p 

//input p,q returns (p-iq)/(p+iq) 
//cos.sin for various angles 
//recursive function 

//tgener(p.twofer) stuffs twofer[2] 
//with the descendants of p 

//returns 1 if both argument edges 
//subtend angles less than SCALE, 
//otherwise returns 0 

//tgener(p.twofer) stuffs twofer 
//with the two descendants of p 

//if both descendants subtend 
//angles less than SCALE, then 
//write two output values and 
//terminate the recursion 

//stuff output values for the first 
//(counter-clockwise) sample point 



ped=twofer, 
pp=-(ped->d); 
qq=ped->c; 
graf[outcount].p=pp; 
graf[outcount].q=qq; 
temp=makecpx(pp,qq); 
ct=temp.x; 
st=temp.y; 

graf[outcount].x=ped->alp.x*ct+ped->beLx*st 

+ped->gam.x+(ped->e.x)*4y(pp*pp+qq*qq); 

graf[outcount].y=ped->aIp.y*ct+ped->bet.y*st 

+ped->gam.y+(ped->e.y)*4y(pp*pp4qq*qq); 

outcount++; //update number of output values 



//stuff output values for the second 
//(counter-clockwise) sample point 



ped++; 
pp=-(ped->d); 
qq=ped->c; 
graf[outcount].p=pp; 
graffoutcount] .q=qq; 
temp=makecpx(pp,qq); 
ct=temp.x; 
st=temp.y; 

graf[outcount].x=ped->alp.x*ct+ped->bet.x*st 

+ped->gam.x; 
graf[outcount].y=ped->alp.y*ct+ped->bet.y*st 

+ped->gam.y; 

outcount++; //update number of output values 



return; 
} 



genertop(twofer,g+ 1); 



//terminate the recursion 



//in the contrary case that one of 
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genertop(twofer+l,§ . .j; 
) 

void generbot(struct edge *p,int g){ 

struct edge twofer[2],*ped; 

struct complex temp, makecpx(int.int); 
double ct,st; 

void generbot(struct edge *,int); 

void bgener(struct edge *,struct edge *); 

int subtendfstruct edge *),pp,qq; 

bgener(p,twofer); 

if(subtend{twofer)) { 



//the descendant edges' subtends ^ 
//large angle, then continue the 
//recursion in counter-clockwise sense 

//recursion for the bottom of 

//the circle, this is entirely 
//analogous to genertop, and 

//only the differences will 
//be noted here 



//differs from genertop in that 
//b,a replaces d,c 



ped=twofer, 
pp=-ped->b; 
qo=ped->a; 
graf[outcountJ.p=pp; 
graffputeount].q=qq; 
temp=makecpx(pp,qq); 
ct=temp.x; 
st=temp.y; 

grafXoutcount].x=-ped->alp.x*ct-ped->bet.x , 'st //differs from 

+ped->gam.x+(ped->e.x)*4y(pp*pp+qq*qq); // genertop in that 

graf[outcount].y=-ped->alp.y*ct-ped->bet.y*st //alp and bet are 

+ped->gam.y+(ped->e.y)*4./{pp*pp+qq*qq); // multiplied by -1 

outcount-H-; 



ped++; 

pp=-ped->b; 

qq=ped->a; 

graf[outcount].p=pp; 

grafloutcount] .q=qq; 

temp=makecpx(pp,qq); 

ct=temp.x; 

st=temp.y; 

graf[outcount].x=-ped->alp.x*ct-ped->bet.x*st 

+ped->gam.x; 
graf[outcount].y=-ped->alp.y*ct-ped->bet.y*st 

+ped->gam.y; 
outcount++; 



//differs from genertop in 
//that b,a replaces d,c 



//differs from genertop 
//in that alp and bet are 
//multiplied by -1 



return; 
} 



generbot(twofer,g+l); 
generbot(twofer+l ,g+l); 



struct complex mult(struct complex u.struct complex v){ 
struct complex temp; 



//multiplication 

//of complex numbers 
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temp.x=u.x*v.x-u.„ 
temp.y-u.x'v.y+v.x* 



struct complex makecpx(int p,int q){ 
struct complex temp; 
double den; 



//input p,q returns 
//the complex number 
//(p-iq)/(p+iq) 



temp.x=(p*p-q*q)/den; 
temp.y=-2*p*q/den; 



int subtend(struct edge *p)( 



nn=(p->a)*(p->c)-Kp->b)*(p->d); 
if(nn*nn<SCAT) 
return(0); 



//subtend returns 1 if each of 
//of the two edges in the 
//array p subtends an angle 
//approximately less than 
//SCALE and returns 0 otherwise 



P++; 

nn=(p->a)*(p->c)+(p->b)*(p->d); 



if<nn*nn<SCAT) 
retum(0); 



void normalin(void){ 

struct complex cz[4],co[4],cm[4]; 



//this routine calculates the modified 
//0.+1.-1 Fourier coeffs czero.cone.cmon 
//from the input Fourier coeffs 



struct complex temp; 



cz{0].x=-l.; 
cz[0].y=0.; 
cz[l].x=0.; 
cz[l].y=0.; 
cz[2].x=-l.; 
cz[2j.y=0.; 
czt3].x=0.; 
cz[3].y=0.; 

co[0].x=0.; 
co[0].y=0.; 
co[l].x=-l.; 
co[l].y=0.; 



//intialize array for czero calculation 



//initialize array for cone calculation 
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//initialize array for cmon calculation 



czer=fourier[0J; 
cone=fourier[l]-, 
cmon=fourier[-l]; 

for(n=2;n<=N;n++){ 

temp=mult(fourier[nJ,cz[n%4]); 

czer.x+=-temp.x; //contribution to czer from the 

czer.y+=-temp.y; . //positive Fourier coeffs 

temp=mult(fourier[-n],cz[(4*N-n)%4]); 



czer.x+=-temp.x; //contribution to czer from the 

czer.y+=-temp.y; //negative Fourier coeffs 

temp=mult(fourier[n],co[n%4]); 

cone.x+=-temp.x; //contribution to cone from the 

cone.y+=-temp.y; //positive Fourier coeffs 



temp=mult(fourier[-n3,co[(4*N-n)%4]); 



cone.x+=-temp.x; 
cone.y+=-temp.y; 



//contribution to cone from the 
//negative Fourier coeffs 



temp=mult(fourier[n],cm[n%4]); 



cmon.x+=-temp.x; 
cmon.y+=-temp.y; 



//contribution to cmon fi 
//positive Fourier coeffs 



from the 



temp=muk(fourier[-n],cm[(4*N-n)%4]); 



cmon.x+=-temp.x; 
cmon.y+=-temp.y; 



//contribution to cmon from the 
//negative Fourier coeffs 



void normalout(void){ 



int n,m; 

struct complex temp,templ,makecpx(int,int); 
double ct,st; 



//this routine alters the 
//output values in the 
//array graf by adding 



//czero+cone*z+cmon*z A (- 1 ), 
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//■where z=(jp-iq)f(p^ . .,; 

for (n=0;n<outcounr,n+-»-){ 

temp=makecpx(graf[n].p,graf[ii3.q); 
templ=mu!t(temp,cone); 

temp.y=-temp.y; 
temp=rnult(temp,cmon); 

graf[n].x+=czer.x+temp 1 ,x+temp.x; 
grafin].y+=czer.y+templ .y+temp.y; 



struct complex cogi(struct edge *p){ //cogi(p) calculates the wavelet coeff for the 

//arrow underlying edge p from the Fourier coeffs 

int a,b,c,d,n; 

struct complex makecpx(int,int), mult(struct complex,struct complex); 
struct complex temp2,templ,temp,crun,cbl,cml,cb2,cm2; 

a=p->a; 

b=p->b; 

c=p->c; 

d=p->d; 

crun.x=0.; 

crun.y=0.; 

cm l=makecpx(b,-a); 
cm2=makecpx(d,-c); 

cbl=cml; 
cb2=cm2; 

for{n=2;n<=N;n++){ 

cbl=mult(cbl,cml); //updated nth power of cml 

cb2=mult(cb2,cm2); //updated nth power of cm2 

if(n<MINMODE)//skip contribution to wavelet coeffs 

continue; //from low-frequency Fourier coeffs 

temp.x=n; 
temp.y=b*d+a*c; 
temp l=mult(cb 1 ,temp); 
temp.y=-temp.y; 
temp2=mul t(cb2, temp); 
temp.x=templ.x+temp2.x; 
temp.y=templ .y+temp2.y; 
templ=mult(fourier[n],temp); 
crun.x+=templ.x; 
crun.y+=templ.y; 

temp.x=-temp.x; 
ternp2=mult(fourier[-n],temp) 
crun.x+=temp2.x; 
crun.y+=temp2.y; 

} 



//contribution to wavelet coeff from 
//nth Fourier coeff for n positive 



//contribution to wavelet coeff from 
//nth Fourier coeff for n negative 
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crun.x=-temp.y/4.; //overall multiple of i/4 

crun.y=temp.x/4. ; 

return(crun); 



void tgener(struct edge *pc,struct edge *pn){ 



//tgener(p,twofer) stuffs twofer 

//with descendants of p for top 
//of circle 



int a,b,c,d,bpd,apc; 



//matrix entries a,b,c,d, 
//bpd=b+d and apc=a+c 



struct complex e,alp,bet,gam; //wavelet coefficient e, trig coefficients alp,bet,gam 



//index 1,2 refers to first.second 
//(counter-clockwise) descendants 
//prefix a,b,c refers to alp,bet,gam 
//update procedure below defines x,y 

//wavelet coeffs of first, second 
//descendant edges 



struct complex axl,ax2; 
struct complex bxl,bx2; 
struct complex cxl,cx2; 
int ayl,ay2,byl,by2,cyl,cy2; 

struct complex eu.et; 

struct complex temp; 



struct complex cogi(struct edge *),mult(struct complex,struct complex); 
struct complex makecpx(int,int); 



double ct,st; 

struct edge *pnp; 

pnp=pn; 
pnp++; 

b=pc->b; 
c=pc->c; 
d=pc->d; 

apc=a+c; 
bpd=b+d; 



//initialize pointers 



//initialize matrix e 



pn->b=b; 

pn->c=apc; 

pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 



x of first descendant 



//stuff matrix of descendant 
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aip=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



//initialize lagged trig coeffs 



e=pc->e; 



//initialize wavelet coeff 



if(a>c){ 



//prepare for stuffing alp,bet,gam in 



a>c 



axl.x=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; 
axl.y=alp.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y; 

ayl=2*(b*b-a*a); 

bxl.x=bet.x+(c*d-a*d-b*c-a*b)*4*e.x; 
bx 1 .y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y; 

byl=4*a*b; 

cx 1 .x=gam.x+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.x; 
cx 1 .y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.y; 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x+4*(d*d-c*c)*e.x; 
ax2.y=alp.y+4*(d*d-c*c)*e.y; 

ay2=2*(c*c-d*d); 

bx2.x=betx+8*e.x*c*d; 
bx2.y=bet.y+8*e.y*c*d; 

by2=-4*c*d; 

cx2.x=gam.x-4*e.x*(c*c+d*d); 
cx2.y=gam.y-4*e.y*(c*c+d*d); 

cy2=2*(c*c+d*d); 
} 

else{ //prepare for stuffing alp,bet,gam in case o=a 

ax 1 .x=alp.x+4*e.x*(a*a-b*b); 
axl.y=aip.y+4*e.y*(a*a-b*b); 

ayl=2*(b*b-a*a); 

bxl.x=bet.x-8*e.x*a*b; 
bxl.y=bet.y-8*e.y*a*b; 

byl=4*a*b; 

cxl.x=gam.x+4*e.x*(a*a+b*b); 
cxl.y=gam.y+4*e.y*(a*a+b*b); 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; 
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ax2.y=al l . , r (a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.y; 
ay2=2*(c*c-d*d); 

bx2.x=bet.x+(a*d-a*b+b*c+c*d)*4*e.x; 
bx2.y=bet.y+(a*d-a*b+b*c+c*d)*4*e.y; 

by2=-4*c*d; 



cx2.x=gam.x+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.x; 
cx2.y=gam.y+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.y; 



cy2=2*(c*c-Ki*d); 
} 



pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 



pn->alp.x=ax 1 .x+ay 1 *et.x; 
pn->alp.y=axl.y+ayl*et.y; 



pn->bet.x=bx 1 .x+by 1 *et.x; 
pn->bet.y =bx 1 .y+by 1 *e t.y ; 



pn->gam.x=cx 1 .x+cy 1 *et.x; 
pn->gam.y=cx 1 .y+cy 1 *et.y ; 



pnp->alp.x=ax2.x+ay2*eu.x; 
pnp->alp.y=ax2.y+ay2*eu.y; 



pnp->bet.x=bx2.x+by2*eu.x; 
pnp->bet.y=bx2.y+by2*eu.y; 



pnp->gam.x=cx2.x+cy2*eu.x; 
pnp->gam.y=cx2.y+cy2*eu.y; 



//calculate wavelet coeffs 
//of the descendants 

//update first a!p 
//update first bet 
//update first gam 
//update second alp 
//update second bet 
//update second gam 



void bgener(struct edge *pc,struct edge *pn){ 



int a,b,c,d,bpd,apc; 

struct complex e,alp,bet,gam; 



//bgener(p,twofer) stuffs twofer 

//with descendants of p for 
//bottom of circle, 
//it is entirely analogous to 
//tgener, and only the 
//differences will be noted here 



struct complex axl,ax2; 
struct complex bxl,bx2; 
struct complex cxl,cx2; 

int ayl,ay2,byl,by2,cyl,cy2; 

struct complex eu,et; 
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struct complex ten. r 



struct complex cogi(struct edge *),mult(struct comple> 
struct complex makecpx(int,int); 



struct edge *pnp; 
double ct,st; 



pnp=pn; 
pnp++; 

a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 



//differs from tgener in that d,-c,-b,a 
//replaces a.b.c.d 



pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 

pnp->a=d; 
pnp->b=-c; 
pnp-xs=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 

e=pc->e; 

if(a>c)[ 



//differs from tgener in that d,-c,-b,a 
//replaces a,b,c,d 



//differs from tgener in that d,-c,-b,a 
//replaces a,b,c,d 



axla=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; 
axl.y=alp.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y; 



bx 1 .x=bet.x+(c*d-a*d-b*c-a*b)*4*e.x; 
bx 1 .y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y ; 



cxljc=gam.x+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.x; 
cxl.y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e,y; 



ax2.x=aip.x+4*(d*d-c*c)*e.x; 
ax2.y=alp.y+4*(d*d-c*c)*e.y; 



ay2=2*(c*c-d*d); 
bx2.x=bet.x+8*e.x*c*d; 
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bx2.y=bei., .-8*e.y*c*d; 
by2=-4*c*d; 

cx2.x=gam.x-4*e.x*(c*c+d*d); 
cx2.y=gam.y-4*e.y*(c*c+d*d); 

cy2=2*(c*c+d*d); 
} 

else{ 

axl .x=alp.x+4*e.x*(a*a-b*b); 
ax 1 .y=alp.y+4*e.y*(a*a-b*b); 

ayl=2*(b*b-a*a); 

bxl.x=bet.x-8*e.x*a*b; 
bxl .y=bet.y-8*e.y*a*b; 

byl=4*a*b; 

cx 1 .x=gam.x+4*e .x* (a *a+b*b); 
cxl .y=gam.y+4*e.y*(a*a+b*b); 

cyl=-2*(a*a+b*b); 

ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; 
ax2.y=alp.y+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.y; 

ay2=2*<c*c-d*d); 

bx2.x=bet.x+(a*d-a*b+b*c+c*d)*4*e.x; 
bx2.y=bet.y+(a*d-a*b+b*c+c*d)*4*e.y; 

by2=-4*c*d; 

cx2.x=gam.x+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.x; 
cx2.y=gam.y+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.y; 

cy2=2*(c*c+d*d); 
} 



pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

pn->alp.x=ax 1 .x+ay 1 *et.x; 
pn->alp.y=ax 1 .y+ay 1 *et.y ; 

pn->bet.x=bx 1 .x+by 1 *et.x; 
pn->bet.y=bx 1 .y+by 1 *et.y ; 

pn->gam.x=cxl .x+cy 1 *et.x; 
pn->gam.y=cx 1 .y+cy 1 *et.y; 

pnp->alp.x=ax2.x+ay2*eu.x; 
pnp->alp.y=ax2.y+ay2*eu .y; 

pnp->bet.x=bx2.x+by2*eu.x; 

68 



WO 00/44104 



PCT/US99/30584 



pnp->bet .y=bx2.y+ ^ „ _* eu .y ; 

pnp->gam.x=cx2.x+cy2*eu.x; 
pnp->gam.y=cx2.y+cy2*eu.y; 

} 



//Here are the subroutines replacing tgener and bgener above in an 
//implementation which takes advantage of renormalization 

void tgener(struct edge *pc,struct edge *pn) 
{ 

int a,b,c,d,bpd,apc; 

struct complex e,alp,bet,gam; 

extern int outcount; 

extern struct pqxy graftCHUNK]; 

struct complex axl,ax2; 
struct complex bxl,bx2; 
struct complex cxl,cx2; 
struct complex eu,et; 
struct complex cogi(struct edge *); 

double ct,st; 
struct edge *pnp; 

pnp=pn; 
pnp++; 
a=pc->a; 
b=pc->b; 
c=pc->c; 
d=pc->d; 

apc=a+c; 
bpd=b+d; 

pn->a=a; 
pn->b=b; 
pn->c=apc; 
pn->d=bpd; 

pnp->a=apc; 
pnp->b=bpd; 
pnp->c=c; 
pnp->d=d; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



e=pc->e; 



69 



WO 00/44104 



PCT/US99/30584 



if(a>c){ 

ax 1 .x=aip.x+4.*e.x; 
ax 1 .y=alp.y+4.*e.y; 
bx I .x=bet.x-4.*e.x; 
bxl .y=bet.y-4.*e.y; 
cxl.x=gam.x; 
cxl.y=gam.y; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2.x=bet.x; 

bx2.y=bet.y; 

cx2.x=gam.x-4.*e.x; 

cx2.y=gam.y-4.*e.y; 

} 



else{ 

ax 1 .x=alp.x+4.*e.x; 

axl.y=alp.y+4.*e.y; 

bxl.x=bet.x; 

bxl.y=bet.y; 

"cx i .x=gam.x+4.*e.x; 

cx 1 .y=gam.y+4.*e.y; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2.x=bet.x+4.*e.x; 

bx2.y=bet.y+4.*e.y; 

cx2.x=gam.x; 

cx2.y=gam.y; 

} 



pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

alp.x=axl.x-2.*et.x; 

bet.x=bxl.x; 

gam.x=cxl.x-2.*et.x; 

alp.y=axl.y-2.*et.y; 
bet.y=bxl.y; 
gam.y=cx 1 .y-2.*et.y; 

pn->aIp.x=(2.*bet.x+alp.x+gara.x)/2.; 

pn->bet.x=(bet.x-alp.x+gam.x); 

pn->gam.x=(2.*bet.x-aIp.x+3.*gam.x)/2.; 

pn->alp.y=(2.*bet.y+alp.y+gam.y)/2.; 

pn->bet.y=(bet.y-alp.y+gam.y); 

pn->gam.y=(2.*bet.y-alp.y+3.*gam.y)/2.; 

aip.x=ax2.x-2.*eu.x; 

bet.x=bx2.x; 

gam.x=cx2.x+2.*eu.x; 

alp.y=ax2.y-2.*eu.y; 
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bet.y=bx2.y; 

ga m .y=cx2.y+2.*eu.y ; 

pnp->alp.x=(-2.*bet.x+alp.x-gam.x)/2.; 
pnp->bet.x=alp.x+bet.x+gam.x; 
pnp->gam.x=(2.*bet.x+alp.x+3 *gam.x)/2.; 

pnp->alp.y=(-2.*bet.y+alp.y-gam.y)/2.; 

pnp->bet.y=alp.y+bet.y+gam.y; 

pnp->gam.y=(2.*bet.y+alp.y+3.*gam.y)/2.; 



double deriv(int a.int b.int c.int d,int p.int q){ 

double x.y.nep; 

x=p*p-q*q; 
y=-2*p*q; 
x=x/(p*p+q*q); 
y=y/(p*p+q*q); 

nep=-(a*a+b*b+c*c+d*d)+x*(a*a+b*b-c*c-d*d)-y*2.*(a*c+b*d); 
retuoi(-2./nep); 



void bgener(struct edge *pc,struct edge *pn) 
i 



int a,b,c,d,bpd,apc; 

struct complex e,alp,bet,gam; 

extern int outcount; 

extern struct pqxy graffCHUNK]; 

struct complex axl,ax2; 

struct complex bxl,bx2; 

struct complex cxl,cx2; 

struct complex eu.et; 

struct complex cogi(struct edge *); 

double ct,sf, 

struct edge *pnp; 

pnp=pn; 
pnp++; 



a=pc->d; 
b=-(pc->c); 
c=-(pc->b); 
d=pc->a; 

apc=a+c; 
bpd=b+d; 

pn->a=bpd; 
pn->b=-apc; 
pn->c=-b; 
pn->d=a; 
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pnp->a=d; 
pnp->b=-c; 
pnp->c=-bpd; 
pnp->d=apc; 

alp=pc->alp; 
bet=pc->bet; 
gam=pc->gam; 



axl.x=alp.x+4.*e.x; 

axl.y=alp.y+4.*e.y; 

bxl.x=bet.x-4.*e.x; 

bxl.y=bet.y-4.*e.y; 

cxl.x=gam.x; 

cxl.y=gam.y; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2.x=bet.x; 

bx2.y=bet.y; 

cx2.x=gam.x-4.*e.x; 

cx2.y=gam.y-4.*e.y; 

} 



else{ 

axl.x=alp.x+4.*e.x; 

ax 1 .y=alp.y+4. *e .y ; 

bxl.x=bet.x; 

bxLy=bet.y; 

cx 1 .x=gam.x+4. *e.x; 

cx 1 .y=gam.y+4. *e.y ; 

ax2.x=axl.x; 

ax2.y=axl.y; 

bx2.x=bet.x+4.*e.x; 

bx2.y=bet.y+4.*e.y; 

cx2.x=gam.x; 

cx2.y=gam.y; 

} 

pn->e=eu=cogi(pn); 
pnp->e=et=cogi(pnp); 

pn->e=eu; 
pnp->e=et; 

alp.x=axl.x-2.*et.x; 

bet.x=bxl.x; 

gam.x=cxl.x-2.*et.x; 

alp.y=ax i .y-2.*et.y; 
bet.y=bxl.y; 
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gam.y=cxl.y-2.*t_. , , 

pn->a!p.x=(2.*bet.x+alp.x+gam.x)/2.; 

pn->bet.x=(bet.x-alp.x+gam.x); 

pn->gam.x=(2.*bet.x-alp.x+3.*gam.x)/2.; 

pn->alp.y=(2.*bet.y+alp.y+gam.y)/2.; 

pn->bet.y=(bet.y-alp.y+gam.y); 

pn->gam.y=(2.*bet.y-a!p.y+3.*gam.y)/2.; 



alp.x=ax2.x-2.*eu.x; 

bet.x=bx2.x; 

gam.x=cx2.x+2.*eu.x; 

alp.y=ax2.y-2.*eu.y; 

bet.y=bx2.y; 

gam.y=cx2.y+2.*eu.y; 

pnp->alp.x=(-2.*bet.x+alp.x-gam.x)/2.; 

pnp->bet.x=alp.x+bet.x+gam.x; 

pnp->garn.x=(2.*bet.x+alp.x+3.*gam.x)/2.; 

pnp->alp.y=(-2.*bet.y+a!p.y-gam.y)/2.; 

pnp->bet.y=alp.y+bet.y+gam.y; 

pnp->gam.y=(2.*bet.y+alp.y+3.*gam.y)/2.; 
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